####################################
# Burglary
####################################

rm(list=ls())

#sink("~/Dropbox/Crime Chile/11_replication/12_RD_figures_burglary_crime.txt")

library(Hmisc)
library(ggplot2)
library(stargazer)
library(foreign)
library(rdrobust)
library(rdd)
library(readstata13)

################
# Prepare data 
################

# reada data
d = read.dta13("~/Dropbox/Crime Chile/11_replication/local_crime_data_chile_2022january.dta")   
names(d)

############################
# RDplots
############################

# Burglary
summary(rdrobust(d$burglary_rate_s,d$margin,cluster = d$cluster,all=TRUE))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA2.pdf", 
          width=7, 
          height=7)
rdplot(y = d$burglary_rate_s, x = d$margin, h= 0.140, nbins = 100, subset = -0.140 <= d$margin & d$margin <= 0.140,
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.541 , 0.238]",
        y.label = "Burglary", x.label = "Margin of victory")
dev.off() 

#################
# Table 4
#################

summary(rdrobust(d$burglary_rate_s,d$margin,cluster = d$cluster, all=TRUE))
pe1 = as.numeric(rdrobust(d$burglary_rate_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv1 = as.numeric(rdrobust(d$burglary_rate_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci1 = as.numeric(rdrobust(d$burglary_rate_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss1 = as.numeric(rdrobust(d$burglary_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$burglary_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess1 = as.numeric(rdrobust(d$burglary_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$burglary_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw1 = as.numeric(rdrobust(d$burglary_rate_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

pe = c(pe1)
pv = c(pv1)
ci_lower = c(ci1[1])
ci_upper = c(ci1[2])
oss = c(oss1)
ess = c(ess1)
bw = c(bw1)
variable = c("Burglary")
data = data.frame(variable,pe,pv,ci_lower,ci_upper,oss,ess,bw)

colnames(data) = c("","Point Estimate","Robust P-value","Robust 95% Confidence Interval lower bound","Robust 95% Confidence Interval upper bound","Overall sample size","Effective sample size","MSE bandwidth")
data

stargazer(data, summary=FALSE, rownames=FALSE, out="~/Dropbox/Crime Chile/11_replication/table4.html")


sink()